home *** CD-ROM | disk | FTP | other *** search
/ TeX 1995 July / TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO / graphics / tiff / tools / tiffmedian.c < prev    next >
C/C++ Source or Header  |  1992-03-26  |  22KB  |  853 lines
  1. #ifndef lint
  2. static char rcsid[] = "$Header: /usr/people/sam/tiff/tools/RCS/tiffmedian.c,v 1.7 92/03/27 14:50:01 sam Exp $";
  3. #endif
  4.  
  5. /*
  6.  * Apply median cut on an image.
  7.  *
  8.  * tiffmedian [-c n] [-f] input output
  9.  *     -c n        - set colortable size.  Default is 256.
  10.  *     -f        - use Floyd-Steinberg dithering.
  11.  *     -lzw        - compress output with LZW
  12.  *     -none        - use no compression on output
  13.  *     -packbits    - use packbits compression on output
  14.  *     -rowsperstrip n    - create output with n rows/strip of data
  15.  * (by default the compression scheme and rows/strip are taken
  16.  *  from the input file)
  17.  *
  18.  * Notes:
  19.  *
  20.  * [1] Floyd-Steinberg dither:
  21.  *  I should point out that the actual fractions we used were, assuming
  22.  *  you are at X, moving left to right:
  23.  *
  24.  *            X     7/16
  25.  *         3/16   5/16  1/16    
  26.  *
  27.  *  Note that the error goes to four neighbors, not three.  I think this
  28.  *  will probably do better (at least for black and white) than the
  29.  *  3/8-3/8-1/4 distribution, at the cost of greater processing.  I have
  30.  *  seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
  31.  *  but I have no idea who the credit really belongs to.
  32.  
  33.  *  Also, I should add that if you do zig-zag scanning (see my immediately
  34.  *  previous message), it is sufficient (but not quite as good) to send
  35.  *  half the error one pixel ahead (e.g. to the right on lines you scan
  36.  *  left to right), and half one pixel straight down.  Again, this is for
  37.  *  black and white;  I've not tried it with color.
  38.  *  -- 
  39.  *                        Lou Steinberg
  40.  *
  41.  * [2] Color Image Quantization for Frame Buffer Display, Paul Heckbert,
  42.  *    Siggraph '82 proceedings, pp. 297-307
  43.  */
  44. #include <stdio.h>
  45. #include "tiffio.h"
  46.  
  47. typedef    unsigned char u_char;
  48. typedef    unsigned short u_short;
  49. typedef    unsigned long u_long;
  50.  
  51. #define    MAX_CMAP_SIZE    256
  52. #define    howmany(x, y)    (((x)+((y)-1))/(y))
  53. #define    streq(a,b)    (strcmp(a,b) == 0)
  54.  
  55. #define    COLOR_DEPTH    8
  56. #define    MAX_COLOR    256
  57.  
  58. #define    B_DEPTH        5        /* # bits/pixel to use */
  59. #define    B_LEN        (1<<B_DEPTH)
  60.  
  61. #define    C_DEPTH        2
  62. #define    C_LEN        (1<<C_DEPTH)    /* # cells/color to use */
  63.  
  64. #define    COLOR_SHIFT    (COLOR_DEPTH-B_DEPTH)
  65.  
  66. typedef    struct colorbox {
  67.     struct    colorbox *next, *prev;
  68.     int    rmin, rmax;
  69.     int    gmin, gmax;
  70.     int    bmin, bmax;
  71.     int    total;
  72. } Colorbox;
  73.  
  74. typedef struct {
  75.     int    num_ents;
  76.     int    entries[MAX_CMAP_SIZE][2];
  77. } C_cell;
  78.  
  79. u_short    rm[MAX_CMAP_SIZE], gm[MAX_CMAP_SIZE], bm[MAX_CMAP_SIZE];
  80. int    bytes_per_pixel;
  81. int    num_colors;
  82. int    histogram[B_LEN][B_LEN][B_LEN];
  83. Colorbox *freeboxes;
  84. Colorbox *usedboxes;
  85. Colorbox *largest_box();
  86. C_cell    **ColorCells;
  87. TIFF    *in, *out;
  88. u_long    rowsperstrip = -1;
  89. u_short    compression = -1;
  90. u_short    bitspersample = 1;
  91. u_short    samplesperpixel;
  92. u_long    imagewidth;
  93. u_long    imagelength;
  94.  
  95. static    get_histogram();
  96. static    splitbox();
  97. static    shrinkbox();
  98. static    map_colortable();
  99. static    quant();
  100. static    quant_fsdither();
  101.  
  102. char    *usage = "usage: tiffmedian [-c n] [-f] [-none] [-lzw] [-packbits] [-rowsperstrip r] input output\n";
  103.  
  104. #define    CopyField(tag, v) \
  105.     if (TIFFGetField(in, tag, &v)) TIFFSetField(out, tag, v)
  106.  
  107. main(argc, argv)
  108.     int argc;
  109.     char **argv;
  110. {
  111.     int i, j, dither = 0;
  112.     char *ifile_name, *ofile_name;
  113.     u_short shortv, config, photometric;
  114.     Colorbox *box_list, *ptr;
  115.     float floatv;
  116.     u_long longv;
  117.  
  118.     num_colors = MAX_CMAP_SIZE;
  119.     for (argc--, argv++; argc > 0 && argv[0][0] == '-'; argc--, argv++) {
  120.         if (streq(argv[0], "-f")) {
  121.             dither = 1;
  122.             continue;
  123.         }
  124.         if (streq(argv[0], "-c")) {
  125.             argc--, argv++;
  126.             if (argc < 1) {
  127.                 fprintf(stderr, "-c: missing colormap size\n%s",
  128.                     usage);
  129.                 exit(1);
  130.             }
  131.             num_colors = atoi(argv[0]);
  132.             if (num_colors > MAX_CMAP_SIZE) {
  133.                 fprintf(stderr,
  134.                    "-c: colormap too big, max %d\n%s",
  135.                    MAX_CMAP_SIZE, usage);
  136.                 exit(1);
  137.             }
  138.             continue;
  139.         }
  140.         if (streq(argv[0], "-none")) {
  141.             compression = COMPRESSION_NONE;
  142.             continue;
  143.         }
  144.         if (streq(argv[0], "-packbits")) {
  145.             compression = COMPRESSION_PACKBITS;
  146.             continue;
  147.         }
  148.         if (streq(argv[0], "-lzw")) {
  149.             compression = COMPRESSION_LZW;
  150.             continue;
  151.         }
  152.         if (streq(argv[0], "-rowsperstrip")) {
  153.             argc--, argv++;
  154.             rowsperstrip = atoi(argv[0]);
  155.             continue;
  156.         }
  157.         fprintf(stderr, "%s: unknown option\n%s", argv[0], usage);
  158.         exit(1);
  159.     }
  160.     if (argc < 2) {
  161.         fprintf(stderr, "%s", usage);
  162.         exit(-1);
  163.     }
  164.     in = TIFFOpen(argv[0], "r");
  165.     if (in == NULL)
  166.         exit(-1);
  167.     TIFFGetField(in, TIFFTAG_IMAGEWIDTH, &imagewidth);
  168.     TIFFGetField(in, TIFFTAG_IMAGELENGTH, &imagelength);
  169.     TIFFGetField(in, TIFFTAG_BITSPERSAMPLE, &bitspersample);
  170.     TIFFGetField(in, TIFFTAG_SAMPLESPERPIXEL, &samplesperpixel);
  171.     if (bitspersample != 8 && bitspersample != 16) {
  172.         fprintf(stderr, "%s: Image must have at least 8-bits/sample\n",
  173.             argv[0]);
  174.         exit(-3);
  175.     }
  176.     if (!TIFFGetField(in, TIFFTAG_PHOTOMETRIC, &photometric) ||
  177.         photometric != PHOTOMETRIC_RGB || samplesperpixel < 3) {
  178.         fprintf(stderr, "%s: Image must have RGB data\n", argv[0]);
  179.         exit(-4);
  180.     }
  181.     TIFFGetField(in, TIFFTAG_PLANARCONFIG, &config);
  182.     if (config != PLANARCONFIG_CONTIG) {
  183.         fprintf(stderr, "%s: Can only handle contiguous data packing\n",
  184.             argv[0]);
  185.         exit(-5);
  186.     }
  187.  
  188.     /*
  189.      * STEP 1:  create empty boxes
  190.      */
  191.     usedboxes = NULL;
  192.     box_list = freeboxes = (Colorbox *)malloc(num_colors*sizeof (Colorbox));
  193.     freeboxes[0].next = &freeboxes[1];
  194.     freeboxes[0].prev = NULL;
  195.     for (i = 1; i < num_colors-1; ++i) {
  196.         freeboxes[i].next = &freeboxes[i+1];
  197.         freeboxes[i].prev = &freeboxes[i-1];
  198.     }
  199.     freeboxes[num_colors-1].next = NULL;
  200.     freeboxes[num_colors-1].prev = &freeboxes[num_colors-2];
  201.  
  202.     /*
  203.      * STEP 2: get histogram, initialize first box
  204.      */
  205.     ptr = freeboxes;
  206.     freeboxes = ptr->next;
  207.     if (freeboxes)
  208.         freeboxes->prev = NULL;
  209.     ptr->next = usedboxes;
  210.     usedboxes = ptr;
  211.     if (ptr->next)
  212.         ptr->next->prev = ptr;
  213.     get_histogram(in, ptr);
  214.  
  215.     /*
  216.      * STEP 3: continually subdivide boxes until no more free
  217.      * boxes remain or until all colors assigned.
  218.      */
  219.     while (freeboxes != NULL) {
  220.         ptr = largest_box();
  221.         if (ptr != NULL)
  222.             splitbox(ptr);
  223.         else
  224.             freeboxes = NULL;
  225.     }
  226.  
  227.     /*
  228.      * STEP 4: assign colors to all boxes
  229.      */
  230.     for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next) {
  231.         rm[i] = ((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2;
  232.         gm[i] = ((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2;
  233.         bm[i] = ((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2;
  234.     }
  235.  
  236.     /* We're done with the boxes now */
  237.     free(box_list);
  238.     freeboxes = usedboxes = NULL;
  239.  
  240.     /*
  241.      * STEP 5: scan histogram and map all values to closest color
  242.      */
  243.     /* 5a: create cell list as described in Heckbert[2] */
  244.     ColorCells = (C_cell **)calloc(C_LEN*C_LEN*C_LEN, sizeof (C_cell *));
  245.     /* 5b: create mapping from truncated pixel space to color
  246.        table entries */
  247.     map_colortable();
  248.  
  249.     /*
  250.      * STEP 6: scan image, match input values to table entries
  251.      */
  252.     out = TIFFOpen(argv[1], "w");
  253.     if (out == NULL)
  254.         exit(-2);
  255.  
  256.     CopyField(TIFFTAG_SUBFILETYPE, longv);
  257.     CopyField(TIFFTAG_IMAGEWIDTH, longv);
  258.     TIFFSetField(out, TIFFTAG_BITSPERSAMPLE, (short)COLOR_DEPTH);
  259.     if (compression != (u_short)-1)
  260.         TIFFSetField(out, TIFFTAG_COMPRESSION, compression);
  261.     else
  262.         CopyField(TIFFTAG_COMPRESSION, compression);
  263.     TIFFSetField(out, TIFFTAG_PHOTOMETRIC, (short)PHOTOMETRIC_PALETTE);
  264.     CopyField(TIFFTAG_ORIENTATION, shortv);
  265.     TIFFSetField(out, TIFFTAG_SAMPLESPERPIXEL, (short)1);
  266.     CopyField(TIFFTAG_PLANARCONFIG, shortv);
  267.     if (rowsperstrip == (u_long)-1)
  268.         rowsperstrip = (8*1024)/TIFFScanlineSize(out);
  269.     TIFFSetField(out, TIFFTAG_ROWSPERSTRIP,
  270.         rowsperstrip == 0 ? 1L : rowsperstrip);
  271.     CopyField(TIFFTAG_MINSAMPLEVALUE, shortv);
  272.     CopyField(TIFFTAG_MAXSAMPLEVALUE, shortv);
  273.     CopyField(TIFFTAG_RESOLUTIONUNIT, shortv);
  274.     CopyField(TIFFTAG_XRESOLUTION, floatv);
  275.     CopyField(TIFFTAG_YRESOLUTION, floatv);
  276.     CopyField(TIFFTAG_XPOSITION, floatv);
  277.     CopyField(TIFFTAG_YPOSITION, floatv);
  278.  
  279.     if (dither)
  280.         quant_fsdither(in, out);
  281.     else
  282.         quant(in, out);
  283.     /*
  284.      * Scale colormap to TIFF-required 16-bit values.
  285.      */
  286. #define    SCALE(x)    (((x)*((1L<<16)-1))/255)
  287.     for (i = 0; i < MAX_CMAP_SIZE; ++i) {
  288.         rm[i] = SCALE(rm[i]);
  289.         gm[i] = SCALE(gm[i]);
  290.         bm[i] = SCALE(bm[i]);
  291.     }
  292.     TIFFSetField(out, TIFFTAG_COLORMAP, rm, gm, bm);
  293.     (void) TIFFClose(out);
  294. }
  295.  
  296. static
  297. get_histogram(in, box)
  298.     TIFF *in;
  299.     register Colorbox *box;
  300. {
  301.     register u_char *inptr;
  302.     register int red, green, blue;
  303.     register u_long j, i;
  304.     u_char *inputline;
  305.  
  306.     inputline = (u_char *)malloc(TIFFScanlineSize(in));
  307.     if (inputline == NULL) {
  308.         fprintf(stderr, "No space for scanline buffer\n");
  309.         exit(-1);
  310.     }
  311.     box->rmin = box->gmin = box->bmin = 999;
  312.     box->rmax = box->gmax = box->bmax = -1;
  313.     box->total = imagewidth * imagelength;
  314.  
  315.     { register int *ptr = &histogram[0][0][0];
  316.       for (i = B_LEN*B_LEN*B_LEN; i-- > 0;)
  317.         *ptr++ = 0;
  318.     }
  319.     for (i = 0; i < imagelength; i++) {
  320.         if (TIFFReadScanline(in, inputline, i, 0) <= 0)
  321.             break;
  322.         inptr = inputline;
  323.         for (j = imagewidth; j-- > 0;) {
  324.             red = *inptr++ >> COLOR_SHIFT;
  325.             green = *inptr++ >> COLOR_SHIFT;
  326.             blue = *inptr++ >> COLOR_SHIFT;
  327.             if (red < box->rmin)
  328.                 box->rmin = red;
  329.                 if (red > box->rmax)
  330.                 box->rmax = red;
  331.                 if (green < box->gmin)
  332.                 box->gmin = green;
  333.                 if (green > box->gmax)
  334.                 box->gmax = green;
  335.                 if (blue < box->bmin)
  336.                 box->bmin = blue;
  337.                 if (blue > box->bmax)
  338.                 box->bmax = blue;
  339.                 histogram[red][green][blue]++;
  340.         }
  341.     }
  342.     free(inputline);
  343. }
  344.  
  345. static Colorbox *
  346. largest_box()
  347. {
  348.     register Colorbox *p, *b;
  349.     register int size;
  350.  
  351.     b = NULL;
  352.     size = -1;
  353.     for (p = usedboxes; p != NULL; p = p->next)
  354.         if ((p->rmax > p->rmin || p->gmax > p->gmin ||
  355.             p->bmax > p->bmin) &&  p->total > size)
  356.                 size = (b = p)->total;
  357.     return (b);
  358. }
  359.  
  360. static
  361. splitbox(ptr)
  362.     register Colorbox *ptr;
  363. {
  364.     int        hist2[B_LEN];
  365.     int        first, last;
  366.     register Colorbox    *new;
  367.     register int    *iptr, *histp;
  368.     register int    i, j;
  369.     register int    ir,ig,ib;
  370.     register int sum, sum1, sum2;
  371.     enum { RED, GREEN, BLUE } axis;
  372.  
  373.     /*
  374.      * See which axis is the largest, do a histogram along that
  375.      * axis.  Split at median point.  Contract both new boxes to
  376.      * fit points and return
  377.      */
  378.     i = ptr->rmax - ptr->rmin;
  379.     if (i >= ptr->gmax - ptr->gmin  && i >= ptr->bmax - ptr->bmin)
  380.         axis = RED;
  381.     else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
  382.         axis = GREEN;
  383.     else
  384.         axis = BLUE;
  385.     /* get histogram along longest axis */
  386.     switch (axis) {
  387.     case RED:
  388.         histp = &hist2[ptr->rmin];
  389.             for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
  390.             *histp = 0;
  391.             for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
  392.                 iptr = &histogram[ir][ig][ptr->bmin];
  393.                 for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  394.                     *histp += *iptr++;
  395.             }
  396.             histp++;
  397.             }
  398.             first = ptr->rmin;
  399.         last = ptr->rmax;
  400.             break;
  401.     case GREEN:
  402.             histp = &hist2[ptr->gmin];
  403.             for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
  404.             *histp = 0;
  405.             for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
  406.                 iptr = &histogram[ir][ig][ptr->bmin];
  407.                 for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  408.                     *histp += *iptr++;
  409.             }
  410.             histp++;
  411.             }
  412.             first = ptr->gmin;
  413.         last = ptr->gmax;
  414.             break;
  415.     case BLUE:
  416.             histp = &hist2[ptr->bmin];
  417.             for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
  418.             *histp = 0;
  419.             for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
  420.                 iptr = &histogram[ir][ptr->gmin][ib];
  421.                 for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
  422.                     *histp += *iptr;
  423.                     iptr += B_LEN;
  424.                 }
  425.             }
  426.             histp++;
  427.             }
  428.             first = ptr->bmin;
  429.         last = ptr->bmax;
  430.             break;
  431.     }
  432.     /* find median point */
  433.     sum2 = ptr->total / 2;
  434.     histp = &hist2[first];
  435.     sum = 0;
  436.     for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
  437.         ;
  438.     if (i == first)
  439.         i++;
  440.  
  441.     /* Create new box, re-allocate points */
  442.     new = freeboxes;
  443.     freeboxes = new->next;
  444.     if (freeboxes)
  445.         freeboxes->prev = NULL;
  446.     if (usedboxes)
  447.         usedboxes->prev = new;
  448.     new->next = usedboxes;
  449.     usedboxes = new;
  450.  
  451.     histp = &hist2[first];
  452.     for (sum1 = 0, j = first; j < i; j++)
  453.         sum1 += *histp++;
  454.     for (sum2 = 0, j = i; j <= last; j++)
  455.         sum2 += *histp++;
  456.     new->total = sum1;
  457.     ptr->total = sum2;
  458.  
  459.     new->rmin = ptr->rmin;
  460.     new->rmax = ptr->rmax;
  461.     new->gmin = ptr->gmin;
  462.     new->gmax = ptr->gmax;
  463.     new->bmin = ptr->bmin;
  464.     new->bmax = ptr->bmax;
  465.     switch (axis) {
  466.     case RED:
  467.         new->rmax = i-1;
  468.             ptr->rmin = i;
  469.             break;
  470.     case GREEN:
  471.             new->gmax = i-1;
  472.             ptr->gmin = i;
  473.             break;
  474.     case BLUE:
  475.             new->bmax = i-1;
  476.             ptr->bmin = i;
  477.             break;
  478.     }
  479.     shrinkbox(new);
  480.     shrinkbox(ptr);
  481. }
  482.  
  483. static
  484. shrinkbox(box)
  485.     register Colorbox *box;
  486. {
  487.     register int *histp, ir, ig, ib;
  488.  
  489.     if (box->rmax > box->rmin) {
  490.         for (ir = box->rmin; ir <= box->rmax; ++ir)
  491.             for (ig = box->gmin; ig <= box->gmax; ++ig) {
  492.                 histp = &histogram[ir][ig][box->bmin];
  493.                     for (ib = box->bmin; ib <= box->bmax; ++ib)
  494.                     if (*histp++ != 0) {
  495.                         box->rmin = ir;
  496.                         goto have_rmin;
  497.                     }
  498.             }
  499.     have_rmin:
  500.         if (box->rmax > box->rmin)
  501.             for (ir = box->rmax; ir >= box->rmin; --ir)
  502.                 for (ig = box->gmin; ig <= box->gmax; ++ig) {
  503.                     histp = &histogram[ir][ig][box->bmin];
  504.                     ib = box->bmin;
  505.                     for (; ib <= box->bmax; ++ib)
  506.                         if (*histp++ != 0) {
  507.                             box->rmax = ir;
  508.                             goto have_rmax;
  509.                         }
  510.                     }
  511.     }
  512. have_rmax:
  513.     if (box->gmax > box->gmin) {
  514.         for (ig = box->gmin; ig <= box->gmax; ++ig)
  515.             for (ir = box->rmin; ir <= box->rmax; ++ir) {
  516.                 histp = &histogram[ir][ig][box->bmin];
  517.                     for (ib = box->bmin; ib <= box->bmax; ++ib)
  518.                 if (*histp++ != 0) {
  519.                     box->gmin = ig;
  520.                     goto have_gmin;
  521.                 }
  522.             }
  523.     have_gmin:
  524.         if (box->gmax > box->gmin)
  525.             for (ig = box->gmax; ig >= box->gmin; --ig)
  526.                 for (ir = box->rmin; ir <= box->rmax; ++ir) {
  527.                     histp = &histogram[ir][ig][box->bmin];
  528.                     ib = box->bmin;
  529.                     for (; ib <= box->bmax; ++ib)
  530.                         if (*histp++ != 0) {
  531.                             box->gmax = ig;
  532.                             goto have_gmax;
  533.                         }
  534.                     }
  535.     }
  536. have_gmax:
  537.     if (box->bmax > box->bmin) {
  538.         for (ib = box->bmin; ib <= box->bmax; ++ib)
  539.             for (ir = box->rmin; ir <= box->rmax; ++ir) {
  540.                 histp = &histogram[ir][box->gmin][ib];
  541.                     for (ig = box->gmin; ig <= box->gmax; ++ig) {
  542.                     if (*histp != 0) {
  543.                         box->bmin = ib;
  544.                         goto have_bmin;
  545.                     }
  546.                     histp += B_LEN;
  547.                     }
  548.                 }
  549.     have_bmin:
  550.         if (box->bmax > box->bmin)
  551.             for (ib = box->bmax; ib >= box->bmin; --ib)
  552.                 for (ir = box->rmin; ir <= box->rmax; ++ir) {
  553.                     histp = &histogram[ir][box->gmin][ib];
  554.                     ig = box->gmin;
  555.                     for (; ig <= box->gmax; ++ig) {
  556.                         if (*histp != 0) {
  557.                             box->bmax = ib;
  558.                             goto have_bmax;
  559.                         }
  560.                         histp += B_LEN;
  561.                     }
  562.                     }
  563.     }
  564. have_bmax:
  565.     ;
  566. }
  567.  
  568. static C_cell *
  569. create_colorcell(red, green, blue)
  570.     int red, green, blue;
  571. {
  572.     register int ir, ig, ib, i;
  573.     register C_cell *ptr;
  574.     int mindist, next_n;
  575.     register int tmp, dist, n;
  576.  
  577.     ir = red >> (COLOR_DEPTH-C_DEPTH);
  578.     ig = green >> (COLOR_DEPTH-C_DEPTH);
  579.     ib = blue >> (COLOR_DEPTH-C_DEPTH);
  580.     ptr = (C_cell *)malloc(sizeof (C_cell));
  581.     *(ColorCells + ir*C_LEN*C_LEN + ig*C_LEN + ib) = ptr;
  582.     ptr->num_ents = 0;
  583.  
  584.     /*
  585.      * Step 1: find all colors inside this cell, while we're at
  586.      *       it, find distance of centermost point to furthest corner
  587.      */
  588.     mindist = 99999999;
  589.     for (i = 0; i < num_colors; ++i) {
  590.         if (rm[i]>>(COLOR_DEPTH-C_DEPTH) != ir  ||
  591.             gm[i]>>(COLOR_DEPTH-C_DEPTH) != ig  ||
  592.             bm[i]>>(COLOR_DEPTH-C_DEPTH) != ib)
  593.             continue;
  594.         ptr->entries[ptr->num_ents][0] = i;
  595.         ptr->entries[ptr->num_ents][1] = 0;
  596.         ++ptr->num_ents;
  597.             tmp = rm[i] - red;
  598.             if (tmp < (MAX_COLOR/C_LEN/2))
  599.             tmp = MAX_COLOR/C_LEN-1 - tmp;
  600.             dist = tmp*tmp;
  601.             tmp = gm[i] - green;
  602.             if (tmp < (MAX_COLOR/C_LEN/2))
  603.             tmp = MAX_COLOR/C_LEN-1 - tmp;
  604.             dist += tmp*tmp;
  605.             tmp = bm[i] - blue;
  606.             if (tmp < (MAX_COLOR/C_LEN/2))
  607.             tmp = MAX_COLOR/C_LEN-1 - tmp;
  608.             dist += tmp*tmp;
  609.             if (dist < mindist)
  610.             mindist = dist;
  611.     }
  612.  
  613.     /*
  614.      * Step 3: find all points within that distance to cell.
  615.      */
  616.     for (i = 0; i < num_colors; ++i) {
  617.         if (rm[i] >> (COLOR_DEPTH-C_DEPTH) == ir  &&
  618.             gm[i] >> (COLOR_DEPTH-C_DEPTH) == ig  &&
  619.             bm[i] >> (COLOR_DEPTH-C_DEPTH) == ib)
  620.             continue;
  621.         dist = 0;
  622.             if ((tmp = red - rm[i]) > 0 ||
  623.             (tmp = rm[i] - (red + MAX_COLOR/C_LEN-1)) > 0 )
  624.             dist += tmp*tmp;
  625.             if ((tmp = green - gm[i]) > 0 ||
  626.             (tmp = gm[i] - (green + MAX_COLOR/C_LEN-1)) > 0 )
  627.             dist += tmp*tmp;
  628.             if ((tmp = blue - bm[i]) > 0 ||
  629.             (tmp = bm[i] - (blue + MAX_COLOR/C_LEN-1)) > 0 )
  630.             dist += tmp*tmp;
  631.             if (dist < mindist) {
  632.             ptr->entries[ptr->num_ents][0] = i;
  633.             ptr->entries[ptr->num_ents][1] = dist;
  634.             ++ptr->num_ents;
  635.             }
  636.     }
  637.  
  638.     /*
  639.      * Sort color cells by distance, use cheap exchange sort
  640.      */
  641.     for (n = ptr->num_ents - 1; n > 0; n = next_n) {
  642.         next_n = 0;
  643.         for (i = 0; i < n; ++i)
  644.             if (ptr->entries[i][1] > ptr->entries[i+1][1]) {
  645.                 tmp = ptr->entries[i][0];
  646.                 ptr->entries[i][0] = ptr->entries[i+1][0];
  647.                 ptr->entries[i+1][0] = tmp;
  648.                 tmp = ptr->entries[i][1];
  649.                 ptr->entries[i][1] = ptr->entries[i+1][1];
  650.                 ptr->entries[i+1][1] = tmp;
  651.                 next_n = i;
  652.                 }
  653.     }
  654.     return (ptr);
  655. }
  656.  
  657. static
  658. map_colortable()
  659. {
  660.     register int *histp = &histogram[0][0][0];
  661.     register C_cell *cell;
  662.     register int j, tmp, d2, dist;
  663.     int ir, ig, ib, i;
  664.  
  665.     for (ir = 0; ir < B_LEN; ++ir)
  666.         for (ig = 0; ig < B_LEN; ++ig)
  667.             for (ib = 0; ib < B_LEN; ++ib, histp++) {
  668.                 if (*histp == 0) {
  669.                     *histp = -1;
  670.                     continue;
  671.                 }
  672.                 cell = *(ColorCells +
  673.                     (((ir>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2) +
  674.                     ((ig>>(B_DEPTH-C_DEPTH)) << C_DEPTH) +
  675.                     (ib>>(B_DEPTH-C_DEPTH))));
  676.                 if (cell == NULL )
  677.                     cell = create_colorcell(
  678.                         ir << COLOR_SHIFT,
  679.                         ig << COLOR_SHIFT,
  680.                         ib << COLOR_SHIFT);
  681.                 dist = 9999999;
  682.                 for (i = 0; i < cell->num_ents &&
  683.                     dist > cell->entries[i][1]; ++i) {
  684.                     j = cell->entries[i][0];
  685.                     d2 = rm[j] - (ir << COLOR_SHIFT);
  686.                     d2 *= d2;
  687.                     tmp = gm[j] - (ig << COLOR_SHIFT);
  688.                     d2 += tmp*tmp;
  689.                     tmp = bm[j] - (ib << COLOR_SHIFT);
  690.                     d2 += tmp*tmp;
  691.                     if (d2 < dist) {
  692.                         dist = d2;
  693.                         *histp = j;
  694.                     }
  695.                 }
  696.             }
  697. }
  698.  
  699. /*
  700.  * straight quantization.  Each pixel is mapped to the colors
  701.  * closest to it.  Color values are rounded to the nearest color
  702.  * table entry.
  703.  */
  704. static
  705. quant(in, out)
  706.     TIFF *in, *out;
  707. {
  708.     u_char    *outline, *inputline;
  709.     register u_char    *outptr, *inptr;
  710.     register u_long i, j;
  711.     register int red, green, blue;
  712.  
  713.     inputline = (u_char *)malloc(TIFFScanlineSize(in));
  714.     outline = (u_char *)malloc(imagewidth);
  715.     for (i = 0; i < imagelength; i++) {
  716.         if (TIFFReadScanline(in, inputline, i, 0) <= 0)
  717.             break;
  718.         inptr = inputline;
  719.         outptr = outline;
  720.         for (j = 0; j < imagewidth; j++) {
  721.             red = *inptr++ >> COLOR_SHIFT;
  722.             green = *inptr++ >> COLOR_SHIFT;
  723.             blue = *inptr++ >> COLOR_SHIFT;
  724.             *outptr++ = histogram[red][green][blue];
  725.         }
  726.         if (TIFFWriteScanline(out, outline, i, 0) < 0)
  727.             break;
  728.     }
  729.     free(inputline);
  730.     free(outline);
  731. }
  732.  
  733. #define    SWAP(type,a,b)    { type p; p = a; a = b; b = p; }
  734.  
  735. #define    GetInputLine(tif, row, bad)                \
  736.     if (TIFFReadScanline(tif, inputline, row, 0) <= 0)    \
  737.         bad;                        \
  738.     inptr = inputline;                    \
  739.     nextptr = nextline;                    \
  740.     for (j = 0; j < imagewidth; ++j) {            \
  741.         *nextptr++ = *inptr++;                \
  742.         *nextptr++ = *inptr++;                \
  743.         *nextptr++ = *inptr++;                \
  744.     }
  745. #define    GetComponent(raw, cshift, c)                \
  746.     cshift = raw;                        \
  747.     if (cshift < 0)                        \
  748.         cshift = 0;                    \
  749.     else if (cshift >= MAX_COLOR)                \
  750.         cshift = MAX_COLOR-1;                \
  751.     c = cshift;                        \
  752.     cshift >>= COLOR_SHIFT;
  753.  
  754. static
  755. quant_fsdither(in, out)
  756.     TIFF *in, *out;
  757. {
  758.     u_char *outline, *inputline, *inptr;
  759.     short *thisline, *nextline;
  760.     register u_char    *outptr;
  761.     register short *thisptr, *nextptr;
  762.     register u_long i, j;
  763.     u_long imax, jmax;
  764.     int lastline, lastpixel;
  765.  
  766.     imax = imagelength - 1;
  767.     jmax = imagewidth - 1;
  768.     inputline = (u_char *)malloc(TIFFScanlineSize(in));
  769.     thisline = (short *)malloc(imagewidth * 3 * sizeof (short));
  770.     nextline = (short *)malloc(imagewidth * 3 * sizeof (short));
  771.     outline = (unsigned char *) malloc(TIFFScanlineSize(out));
  772.  
  773.     GetInputLine(in, 0, goto bad);        /* get first line */
  774.     for (i = 1; i < imagelength; ++i) {
  775.         SWAP(short *, thisline, nextline);
  776.         lastline = (i == imax);
  777.         GetInputLine(in, i, break);
  778.         thisptr = thisline;
  779.         nextptr = nextline;
  780.         outptr = outline;
  781.         for (j = 0; j < imagewidth; ++j) {
  782.             int red, green, blue;
  783.             register int oval, r2, g2, b2;
  784.  
  785.             lastpixel = (j == jmax);
  786.             GetComponent(*thisptr++, r2, red);
  787.             GetComponent(*thisptr++, g2, green);
  788.             GetComponent(*thisptr++, b2, blue);
  789.             oval = histogram[r2][g2][b2];
  790.             if (oval == -1) {
  791.                 int ci;
  792.                 register int cj, tmp, d2, dist;
  793.                 register C_cell    *cell;
  794.  
  795.                 cell = *(ColorCells +
  796.                     (((r2>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2) +
  797.                     ((g2>>(B_DEPTH-C_DEPTH)) << C_DEPTH ) +
  798.                     (b2>>(B_DEPTH-C_DEPTH))));
  799.                 if (cell == NULL)
  800.                     cell = create_colorcell(red,
  801.                         green, blue);
  802.                 dist = 9999999;
  803.                 for (ci = 0; ci < cell->num_ents && dist > cell->entries[ci][1]; ++ci) {
  804.                     cj = cell->entries[ci][0];
  805.                     d2 = (rm[cj] >> COLOR_SHIFT) - r2;
  806.                     d2 *= d2;
  807.                     tmp = (gm[cj] >> COLOR_SHIFT) - g2;
  808.                     d2 += tmp*tmp;
  809.                     tmp = (bm[cj] >> COLOR_SHIFT) - b2;
  810.                     d2 += tmp*tmp;
  811.                     if (d2 < dist) {
  812.                         dist = d2;
  813.                         oval = cj;
  814.                     }
  815.                 }
  816.                 histogram[r2][g2][b2] = oval;
  817.             }
  818.             *outptr++ = oval;
  819.             red -= rm[oval];
  820.             green -= gm[oval];
  821.             blue -= bm[oval];
  822.             if (!lastpixel) {
  823.                 thisptr[0] += blue * 7 / 16;
  824.                 thisptr[1] += green * 7 / 16;
  825.                 thisptr[2] += red * 7 / 16;
  826.             }
  827.             if (!lastline) {
  828.                 if (j != 0) {
  829.                     nextptr[-3] += blue * 3 / 16;
  830.                     nextptr[-2] += green * 3 / 16;
  831.                     nextptr[-1] += red * 3 / 16;
  832.                 }
  833.                 nextptr[0] += blue * 5 / 16;
  834.                 nextptr[1] += green * 5 / 16;
  835.                 nextptr[2] += red * 5 / 16;
  836.                 if (!lastpixel) {
  837.                     nextptr[3] += blue / 16;
  838.                         nextptr[4] += green / 16;
  839.                         nextptr[5] += red / 16;
  840.                 }
  841.                 nextptr += 3;
  842.             }
  843.         }
  844.         if (TIFFWriteScanline(out, outline, i-1, 0) < 0)
  845.             break;
  846.     }
  847. bad:
  848.     free(inputline);
  849.     free(thisline);
  850.     free(nextline);
  851.     free(outline);
  852. }
  853.